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1 .  INTRODUCTION 


An  important  problem  in  turbomachinery  is  the  prediction  of 
the  flutter  boundaries  of  the  compressor.  In  order  to  compute 
these  boundaries,  the  unsteady  aerodynamic  forces  need  to  be 
understood.  The  aerodynamic  causes  of  flutter  in  a  compressor 
can  be  very  complex,  for  examp]  ;,  the  interaction  between  the 
flows  induced  by  the  rotor  and  stator.  However,  an  important 
class  of  flutter  is  caused  by  a  relatively  simple  aerodynamic 
flow  due  to  the  vibration  of  the  blades  on  a  given  compressor  row 
without  the  added  complexity  of  exterior  interactions.  This 
latter  problem  is  the  simplest  unsteady  aerodynamic  phenomenon  of 
interest.  In  a  compressor  the  oncoming  velocity  varies  along  the 
compressor  blade,  increasing  from  subsonic  speeds  at  the  hub  to 
perhaps  transonic  speeds  at  the  tip;  it  is  the  prediction  of 
unsteady  transonic  flow  in  a  compressor  that  is  the  most 
difficult.  This  report  is  concerned  with  the  prediction  of 
flutter  in  a  vibrating  compressor  row  when  the  flow  is  transonic. 

The  work  reported  here  is  concerned  with  solving  the  classic 
high  frequency  unsteady  transonic  small  disturbance  equation  for 
cascade  flows.  A  second  topic  is  the  validation  of  the  cascade 
theory  of  Reference  1. 

A  major  problem  with  potential  theory  is  that  it  is  not 
valid  (Reference  2)  when  shock  waves  are  present  and,  as  a 
consequence,  shock  locations  are  not  predicted  adequately  for 
medium  to  strong  shock  strengths,  although  the  agreement  for 
flows  with  weak  shocks  is  acceptable.  This  disadvantage  of  the 
potential  formulation  can  be  overcome  by  modifying  the  theory  to 
give  the  correct  shock  jump.  This  has  been  done  by  Nixon 
(Reference  3)  for  the  steady  transonic  small  disturbance  (TSD) 
equation  and  by  Kerlick,  Nixon,  and  Ballhaus  (Reference  4)  for 
the  unsteady  low  frequency  equation.  An  extension  of  the 
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work  in  Reference  4  is  reported  in  Reference  5  which  is  attached 
as  Appendix  1. 

Most  of  the  numerical  methods  for  predicting  unsteady 
transonic  flow  stem  from  the  work  of  Ballhaus  and  Goorjian 
(Reference  6),  who  developed  an  implicit  algorithm  to  solve  the 
low  frequency  TSD  equation  for  isolated  airfoils.  An  extension 
of  this  algorithm  to  high  frequency  motion  was  made  by  Rizzetta 
and  Chin  (Reference  7).  For  cascade  flows  a  version  of  the 
Rizzetta-Chin  algorithm  was  used  by  Kerlick  and  Nixon 
(Reference  8)  to  develop  a  method  for  cascade  flows  for 
unstaggered  cascades.  An  alternative  to  the  use  of  the  unsteady 
TSD  equation  has  been  considered  by  Verdon  and  Caspar 
(Reference  9)  who  use  a  time  linearized  perturbation  of  the  full 
potential  equation  for  staggered  and  unstaggered  cascades  in  a 
subsonic  flow. 

In  a  cascade  flutter  analysis  the  dominant  parameter  is  the 
interblade  phase  angle  which  can  determine  the  conditions  at 
which  the  blade  will  flutter.  In  a  numerical  method,  such  as 
those  discussed  earlier,  a  flutter  calculation  can  involve 
computing  a  test  case  for  each  value  of  the  interblade  phase 
angle.  This  is  computationally  expensive  and,  in  order  to  reduce 
the  cost,  a  technique  developed  by  Nixon  (Reference  1)  may  be 
used.  This  technique,  which  is  based  on  indicial  theory, 
decouples  the  cascade  problem  into  a  series  of  elementary 
problems,  which  are  parametrically  related,  and  only  one  solution 
need  be  computed.  This  elementary  solution  consists  of  a  number 
of  blades,  all  stationary  with  the  exception  of  one  blade  which 
undergoes  a  step  change.  Once  computed,  the  full  unsteady 
solution  for  any  frequency  and  interblade  phase  angle  can  be 
constructed  by  a  judicious  superposition. 


In  a  previous  report  (Reference  10)  the  results  obtained  by 
the  cascade  code  described  in  Reference  8  and  the  indicial  method 
described  in  Reference  1  are  examined.  Several  areas  of  concern 
arise  in  this  examination.  First,  for  a  flat  plate  at  subsonic 
conditions,  the  method  of  Kerlick  and  Nixon  Reference  8  did  not 
agree  very  well  with  the  results  of  Verdon  and  Casper  (Reference 
9).  Second,  the  indicial  method  did  not  give  the  super-resonance 
phenomena  noted  by  Verdon  and  Casper.  It  is  desirable  that  these 
points  be  clarified  before  further  work  is  done.  The  question  of 
super-resonance  on  the  indicial  method  has  been  answered  by 
including  a  (nominally)  infinite  number  of  blades  in  the 
cascade.  This  is  done  analytically  rather  than  computationally. 

In  the  above  studies  some  questions  regarding  the  far  field 
boundary  conditions  used  in  cascade  computations  are  raised.  It 
is  possible  that  the  differences  in  the  boundary  conditions 
account  for  the  difference  in  the  flat  plate  results  noted 
above.  A  discussion  is  given  in  the  text. 

In  order  to  compute  flows  with  strong  shock  waves  the 
techniques  of  References  4  and  5  have  been  incorporated  into  the 
computer  code.  The  details  are  given  in  Appendix  1. 

The  final  issue  to  be  addressed  in  the  report  is  that  of 
nonunigue  solutions.  It  has  been  known  for  some  years  that  the 
potential  equation  gives  nonunique  solutions  (References  11,  12) 
for  isolated  airfoils.  It  is  established  here  that  these 
solutions  also  occur  for  steady  cascade  flows.  For  unsteady 
flows  the  unsteady  perturbation  seems  to  be  unique  although  the 
perturbation  is  around  an  artificial  value  of  lift.  These 
computations  are  for  a  reduced  frequency  of  unity  and  is  in 


agreement  with  the  results  of  Williams  et  al.  (Reference  13)  for 
isolated  airfoils.  However,  for  low  frequencies  an  artificial 
behavior  is  evident  for  isolated  airfoils  (Reference  13)  and  may 
occur  for  cascade  flows. 

The  analyses  of  Steinhoff  and  Jameson  (Reference  11),  Salas 
et  al.  (Reference  12),  and  Williams  et  al.  (Reference  13)  are 
essentially  numerical  experiments.  Nixon  (Reference  14)  has 
attempted  a  more  analytic  study  of  the  problem.  Reference  14  is 
atcached  as  Appendix  2.  In  Reference  14  Nixon  has  shown  that  for 
isolated  airfoils  there  is  the  possibility  of  multiple  solutions 
due  to  the  presence  of  eigensolutions  in  the  eguations.  There 
does  not  seem  to  be  any  way  to  avoid  these  solutions;  in  fact  it 
is  indicated  that  a  stable  algorithm  helps  the  appearance  of  the 
eigensolutions.  In  the  present  work  the  analysis  of  Reference  14 
is  extended  to  include  steady  cascade  effects.  It  is  found  that 
the  effect  of  the  eigensolutions  varies  as  the  inverse  square  of 
the  gap.  It  is  likely  that  for  practical  cascade  flows  the 
nonunique  solutions  do  not  arise. 

2.  COMMENTS  ON  THE  ACCURACY  OF  THE  NEAR  CASCADE  CODE 

In  the  previous  work  on  this  contract,  reported  in  Reference 
10  a  discrepancy  between  the  results  of  the  computer  code  CTRAN, 
developed  under  Contract  N00018-81-C-0169 ,  and  the  code  of  Verdon 
and  Casper  (Reference  9)  for  the  subsonic  flow  through  a  cascade 
of  flat  plates  was  noted.  Since  the  flat  plate  is  a  classical 
case  this  is  a  disturbing  error.  Hence  some  thought  has  been 
given  as  to  the  reasons  for  the  disagreement. 

For  subsonic  flows  the  Verdon  and  Caspar  equation  is 
equivalent  to  the  small  disturbance  equation 
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where  0  is  the  angle  the  oncoming  stream  makes  with  the  x-axis. 
The  equation  used  by  Nixon  and  Kerlick  (References  8  and  10  )  is 
the  classic,  isolated  airfoil  equation 
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Hence,  it  can  be  seen  that  the  Verdon  and  Caspar  equations 
include  the  extra  terms  due  to  the  y  derivatives  which  are 
necessary  to  model  highly  curved  blades.  The  boundary  conditions 
for  both  equations  are  that  there  is  no  flow  through  a  solid 
boundary.  For  the  case  of  a  flat  plate  at  angle  of  attack  a  the 
boundary  condition  is 

♦y(x,±  0)  ■-<*  +  £-  (x-xQ)  afc  (3) 

where  xQ  is  some  pitching  axis. 

If  0  is  small,  Equations  (1)  and  (2)  become  identical.  In 
the  comparisons  in  Reference  10  0  is  zero  and,  thus,  the 
discrepancy  reported  cannot  be  due  to  the  different  equations.  A 
possible  cause  of  the  discrepancy  may  be  the  different  far  field 
boundary  conditions  (see  Section  3)  but  this  has  not  been 
investigated. 

3.  FAR  FIELD  BOUNDARY  CONDITIONS 


The  far  field  boundary  condition  used  by  Verdon  and  Caspar 
is  that  t  has  the  asymptotic  form  upstream  of  the  cascade 
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(4) 


where 

q.  =  (2nj+o)/h,  6.=  M  ( w  +  U  sin  8  g . ) 

3  3  •  “  3 

O2  -  a^Jl  "  cos20),  d2  =  (02q.2-  62)M 2a<o2  (5) 

and  a  is  the  speed  of  sound  in  the  undisturbed  flow,  a  is  the 
interblade  phase  angle  and  u  is  the  frequency  of  oscillation. 

x_w  denotes  the  upstream  computational  boundary  and  h  is  the 
cascade  gap. 

The  constant  bj  is  given  by 

bj  =  h  1/)J+h  4>(  5/n)exp(-qjn)dn  (6) 

Far  upstream  of  the  cascade  a  specified  boundary  condition  is 
that  the  velocity  components  are  given. 

Thus  as  x  ♦  -« 


<)>(x,n)  n  sin  8 


(7) 


In  the  limit  as 


and,  from  equations  (6)  and  (7) 


b .  ♦  U  sin  0 
3 


In  the  study  of  asymptotic  behavior  there  are  two  limits  that  can 
be  examined  as  the  gap  approaches  infinity.  First,  if  the  point 
x  approaches  -®  faster  than  h  ♦  ®  then  the  effect  at  x_w  will 
always  be  that  of  a  cascade.  Second,  if  the  point  x  does  not 
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approach  infinity  faster  than  h  +  «  then  the  behavior  at  x  will 

- •  00 

be  that  of  an  isolated  airfoil.  Mathematically  the  first  case  is 
when 

&i?  s-  -tt  *  °  <»' 

while  the  second  case  is  when 


(10) 


In  the  case  of  an  isolated  airfoil  the  classic  limit,  on 
y  =  0  as  x  ♦  -•  ,  is  (see,  for  example,  Reference  15) 
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As  x  ♦  -«  the  value  of  <t>(x,y)  given  by  Eauation  (11)  decays  as 
| x | /2  whereas  the  form  given  by  Equation  (4)  behaves  as 

1  o 

I  U  sin  0  exp  (i[M^n  ±  K)  (x-x  )} 

W  00  ““GO 


and  b.  -*■  constant  as  h  +  «.  Note  that  as  h  +  •  the  number  of 
3 

blade  gaps  must  reduce  to  two  with  j  =  ±  1  denoting  the  far  field 
blades  at  n  =  ±  “• 


Equation  (13)  does  not  give  the  isolated  airfoil  result 
as  h  +  «  if  the  limit  of  Equation  (10)  is  imposed.  It  is 
suggested  that  this  limit  is  the  correct  one  to  use  for  a  far 
field  computational  boundary  since  x_w  must  be  finite  to  be  of 
any  use. 
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STRONG  SHOCK  THEORY 


In  Reference  5  a  modification  to  the  unsteady  high  frequency 
small  disturbance  equation  is  described  which  alters  the  shock 
strength  to  give  the  Rankine  Hugoniot  value.  These  types  of 
modifications  improve  the  accuracy  of  the  small  disturbance 
equation.  Reference  5  is  attached  as  Appendix  1. 

The  theory  of  Reference  5  has  been  applied  to  the  cascade 
code  CTRAN.  Since  CTRAN  and  the  base  code  XTRAN2.L  used  in 
Reference  5  are  similar  the  mechanics  of  coding  the  appropriate 
modifications  proved  to  be  easy.  A  detailed  description  of  the 
modifications  is  given  in  Appendix  1.  Typical  results  for  a 
cascade  are  shown  in  Figures  1  and  2.  The  "strong  shock" 
correction  for  a  cascade  is  similar  to  that  for  isolated  airfoil 
in  that  the  shock  waves  are  weakened  and  move  forward.  The 
section  is  a  NACA0012;  the  amplitude  is  1/4°,  M<jo=  0.845  and  h  = 
30.  The  large  gap  is  used  to  indicate  the  effect  of  the 
modification.  For  more  realistic  values  of  the  gap  the  strong 
shock  modification  may  prove  to  be  unnecessary  because  the  shock 
waves  in  such  conditions  are  relatively  weak  anyway  unless  the 
flow  is  nearly  choked. 

5.  COMMENTS  ON  THE  INDICIAL  THEORY  FOR  CASCADE  FLOWS 


In  Reference  1  a  method  of  computing  the  unsteady  flow 
through  a  cascade  for  any  interblade  phase  angle  using  only  one 
set  of  computational  data  is  described.  The  computational  data 
are  calculated  for  a  cascade  with  only  one  blade  moving  and  the 
others  held  stationary.  A  further  simplification  is  that  only 
the  indicial  response  for  the  cascade  needs  be  computed;  from 
this  result  the  flow  due  to  arbitrary  time  dependent  motions  can 
be  constructed.  In  the  previous  report  on  the  contract 
(Reference  10)  attempts  to  compute  the  unsteady  flow  through 
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where  CL  is  the  value  of  CL  as  h  +  «,  hQ  is  some  value  of  h  to 
be  determined  and  A  is  a  constant. 

The  computer  code  XTRAN2L,  used  for  isolated  airfoil 
problems,  was  modified  to  treat  a  steady  cascade  by  the  addition 
of  periodic  boundary  conditions.  A  series  of  calculations  for  a 
NACA0012  cascade  at  Mb  =  0.845  and  a  =  0  was  performed  with 
various  values  of  the  gap,  h.  The  results  are  shown  in  Figure  6. 

If  the  value  is  h  when  CL  is  zero  (h  =  22)  and  the  value  of 
CL  at  h  *  25  are  used  to  compute  hQ  and  A  in  Equation  (60)  then 
CL  can  be  computed.  This  value  is  also  shown  in  Figure  4  and 
agrees  very  well  with  the  computed  results.  This  good  agreement 
indicates  that  the  compatibility  equation.  Equation  (55), 
represents  the  behavior  of  nonunique  solutions  and  that  the 
general  conclusions  given  in  Appendix  2  apply  to  cascade  flows. 

6.  CONCLUDING  REMARKS 

A  number  of  points  arising  from  previous  work  on  cascade 
flows  have  been  investigated  and  most  have  been  clarified.  There 
is  still  some  question  as  to  the  correct  far  field  boundary 
conditions  to  apply  to  unsteady  cascade  computations.  However, 
the  existence  of  nonunique  solutions  to  the  potential  flow  around 
cascade  has  been  verified  and  examined. 
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If  the  gap  is  large  the  computationally  effective  flow  field 
will  be  similar  to  that  of  an  isolated  airfoil  since  u  and  Au 
will  approach  zero  as  n  ♦  °°.  Consequently,  the  main 
contribution  to  the  field  integral  will  come  only  from  values  of 
u  up  to  some  distance  h  from  the  blades.  In  other  words 

■/JJ- f<n)dn  ~  j£f(n)dn  (57) 

where  h  is  independent  of  h.  Consequently  the  field  integrals  in 
Equations  (55)  and  (52)  can  be  written,  to  a  first  approximation. 


l  LI  (5af-)2ids1  =  /:  /soti(9l«f 


n=o  1 


1  *1’  *f, 


+  i  H  (x  ,y ) 
n  1 


Z  ~  i  T.  U  I  Ak| 


fldSl  +  ~2  H2(x'y) 
h 


where  (x,y)  and  H2  (x,y)  are  independent  of  h  and  f^  is 
similar  to  the  eigensolution  for  the  isolated  airfoil,  g-i  and 

A  A 

If^  are  associated  with  fj.  is  a  subset  of  Sj  with  a 

range  0  <  n  <  h  rather  than  0  <  <  h. 

By  examining  Equations  (55),  (58)  and  (59)  it  can  be  deduced 

that  the  compatibility  equation.  Equation  (55),  is  similar  to 

that  of  the  isolated  airfoil  case,  with  the  addition  of  a 
perturbation  depending  on  the  parameter  — From  the  general 
transonic  perturbation  theory  of  Reference  16  it  is  expected  that 
the  variation  of  CL  in  a  nonunique  solution  will  vary  as  — 


If  the  perturbation  parameter  is  of  the  form  — 2  it  follows 
that  the  lift  coefficient  for  a  cascade  is  of  the  £orm 
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(X) 


Consider  now  the  case  when  the  gap  is  large. 

To  a  first  approximation 

uT  ( x,y )  *  u,J,o)(x,y)  +  ~2  u^tx.y)  (56) 

L  L  h  L 

where  u*°^  is  the  value  for  the  isolated  airfoil  and  u^1 *  is  a 
1 L  *  r 

function  independent  of  h. 
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f[o)  -  fn  (5,0) 
n  n 


An  alternate  form  of  Equation  (42)  is 


F  +  F  *  T  f 
xx  yy  nto  nx 

with  the  boundary  condition 


Fy(x ,±nh)  -  0, 


where 


F(x,y)  «  Jxr~ —  {f.(C,y) 

u(c,y)  1 


I1  ~u VI  f(o)(5)}de 

t 1 — u ( 5,o)]  l 


*  /x  { Au  (  £ ,  y )  -  Au(C,o)}dC  (4( 

Equation  (44)  and  (46)  are  similar  to  Equation  (44)  and  (46)  in 
Reference  14  and  the  nonuniqueness  analysis  is  then  similar  to 
that  in  Reference  14. 

The  symmetric  part  of  Equation  (37)  is 


“  -  ■  “tl  -  j„ 


where 


Ku’  -  i  "'u’lx-y)  + 


*  i  i  Av  (x  -  5) 

uT  (x,y)  *  l  2ir  ^o  - 2 - T 

TL  n  =  —  2ir  °  t(x-E)2+  (y-nh)2] 
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The  asymmetric  part  of  Equation  (37)  can  be  written  as 


“  ,  Av  (x-C) 

Au  -Auu  =  I  J  /  t - j - 2“ 

n=-«  [{x-£)^+  (y-nh)^) 


1  ^  AuQ(y+nh) 

C  2^°  t(x-C)2+(y+nh)2]' 


~  l  lo  -  Kp”lx,-y)J  Auuds, 

n*0  bl  5X  5X  1 


From  Equation  (28)  and  (30) 


Auu-(2n+l)Au  «  f  (39) 

o  n 

This  equation  specifies  Au  in  the  flow  field  in  terms  of  fn. 

Since  Equation  (38)  also  specifies  Au  in  the  flow  field  in  terms 
of  Auq  there  must  be  a  compatibility  relation. 

Equation  (38)  can  be  written  as 

Au  -  Auu  *  -  l  t  /c  /  AK  *n) [Auu-(2n+l)Au  IdS.  (40) 

n*0  **  &1  tx  o  i 


where 


tK^^Xry)  -  K^n)(x,-y)] 


Using  Equation  (39),  Equation  (40)  can  be  written  as 


Z  i±:~-  f{0)>  -  ff  l  jJs  ndSl 

u  1  (l-uQ)  1  1  n-04  S1  n 


°  Uo-  1  f(o) 

Au  -  I  {f .+  ] 

u  (u0-  1) 
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where 


L(vq)  is  an  integral  operator  on  vq 


such  that 


1 

2* 


00 

l 


L(vq)(x-5) 

t(x-e)2+  n2h2] 


<35 


v 

o 


(33) 


The  fact  that  a  lifting  solution  to  the  linear  problem  exists 
such  that 


Au  «  L (v  )  (34) 

o  o 

is  significant  evidence  that  L(vq)  exists. 

Equation  (32)  can  be  written  in  the  form 

0  -  -  l  /_  {Au(C,n)uU,n)  -  ( 2u+l ) [ Au  -L( v) ] }dS_  (35) 

n=0  &1  4yo  °  1 

Equation  (44)  is  similar  to  Equation  (30)  and  indicates  that  a 
nonzero  function,  fn,  exists,  which  satisfies  Equation  (30),  and 
is  given  by 


f  *  au( 5,n)u( c, n)-( 2u+l ) [ Au  -  L ( v  ) ] 
n  o  o 


(36) 


where  the  right-hand  side  is  generated  from  a  real  transonic 
solution. 


Now  differentiate  Equation  (21)  with  respect  to  x,  giving 
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U  r  ,1 

u"  2  "  L  {2*'o 


•1  r" 


Avo(x-C) 


((x-S)2+  (y-nh)2) 


<35  +  L  ! 


AuQ(y-nh) 


2li  °  [  (x-£)  2+  (y-nh)2] 


<30 


- 1  isi*l*u2'ss 


(37) 
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where 


v(x,o)  =  ^  [v(x,+o)  +  v(x-0 ) 1 


and  it  has  been  noted  that,  because  of  periodicity. 


Au(?,nh)  =  Au(£,o)  =  Au 


The  subscript  o  denotes  a  value  on  y  =  0.  Equation  (25)  can  be 
written  as 


v(x,o)  =  -  l  Jc  fKr„n  tAu(c,n)u(C,n) 
n=o  &1  *yo 


(2n+l)Au  ]dS, 
o  1 


where  now  covers  one  blade  gap 

K(n)=  •iirln((x-C)2+  (n-nh-y)2]  (29) 

Equation  (28)  is  similar  to  one  in  Reference  14  and  for  a  lifting 
solution  to  exist  when 


v(x,o)  =  0 


a  function,  fn,  must  exist  such  that 


£n 


In  order  to  establish  the  existence  of  fn  it  is  necessary  to 
write  Equation  (25)  as 


»  -  J  <5,  / 


( Au  -L  ( v  ) ]  (x-5)  ,  2 

o  - 2 - 1 - TT~  de}  "  lJsfKcy  u  ds 

t(x-e)  +  nzhz)  2v  S  *yo 
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Analysis  of  Nonunique  Solutions 

A  cascade  is  composed  of  an  infinite  number  of  blades  or 
airfoils  and  an  integral  representation  of  the  TSD  equation  can 
be  constructed  in  a  similar  way  as  that  for  an  isolated  airfoil 
(described  in  Reference  14)  using  the  principle  of 
superposition.  Thus, 

ao 

♦  (x,y)  =  l  {^"[Avic£n)-  A^"5  )d5}-l/2/s/K  u2dS  (21) 

n=— •  o 

where 


K^n>  -  jutnj  (x-02+  (nh-y)2j 


(22) 


and 


K  =  jn*n| (x-5)2+  (y-n) 2 |  (23) 

Differentiation  of  Eauation  (21)  with  respect  to  y  gives 

v(x.y)  -  £  (J-  ^U,nh)(y-nh) 

n=—  [(x-£)^+  (y-nh)  ] 

-  T  W*48 

The  operator  "A"  is  defined  as 

Af(x,y)  =  \  [f (x,+y)  -  f(x-y))  .  (24) 


-  r  .  ^(^.,nh)  (x-o  dc} 

[  (x-£) (y-nh) z] 


On  y  s  0  the  first  sum  yields  the  value  note  that  since 

Av(5,n)  is  periodic  about  each  blade  the  terms  for  n  *  ±  N 
cancel.  Equation  (24)  can  thus  be  written  as 


v(x,o) 


-  I 


n  =  -«* 


Au  (  C  , o )  ( x-Q 
2  2  2 
[(x-C)  +  nTO 


<35 


1  W/45 


(25) 


addition  to  the  cascade  result,  results  for  both  a  free  air  and 
for  a  solid  wall  boundary  condition  are  shown;  in  these  cases  the 
far  field  boundary  is  at  the  same  location  as  for  the  cascade 
result.  It  can  be  seen  that  the  behavior  with  the  various 
boundary  conditions  is  different  even  though  the  boundary  is  at  a 
nominally  far  field  location.  The  behavior  of  a  nonunique 
solution  with  decreasing  gap  is  shown  in  Figure  6  and  is  compared 
with  the  result  for  a  solid  wall  boundary  condition.  The 
difference  is  quite  significant.  It  also  appears  that  the 
cascade  effect  inhibits  the  appearance  of  the  nonuniqe  solution 
since  it  does  not  appear  for  practical  gap/chord  ratios  for  this 
particular  section. 

The  cascade  code  was  used  to  investigate  the  behavior  of  the 
nonunigue  solutions  in  oscillatory  flows.  These  oscillatory 
solutions  were  started  from  either  a  zero  lift  condition  or  a 
nonunique  steady  state.  The  results  are  shown  in  Figure  7;  and 
it  can  be  seen  that  the  oscillatory  perturbation  is  nearly 
sinusoidal  but  is  about  a  non  zero  lift.  The  section  is  the 
NACA0012,  *  0.845,  the  pitch  amplitude  is  1/4°  and  the  reduced 

frequency  is  1.0.  It  is  also  noted  that  the  solution  starting 
from  a  zero  lift  condition  is  not  converged  but  appears  to  be 
approaching  the  other  results.  As  a  final  result  the  cascade 
code  with  the  strong  shock  conditions  was  run  for  the  same  case 
and  the  results  seem  to  be  of  the  unique  variety.  The  steady 
state  solution  is  symmetric  in  this  case  in  order  to  provide 
symmetric  controls  for  the  strong  shock  addition;  asymmetry  is 
introduced  by  allowing  a  non  zero  mean  angle  of  attack  for  the 
first  200  time  steps. 

It  is  desirable  to  examine  the  behavior  of  these  nonunique 
solutions  and  this  is  done  in  the  next  section. 
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CL  will  be  infinite.  For  the  flat  plate  reported  case  by  Verdon 
and  Caspar  with  «  0. 5,  h  =  1 ,  v  =  1  Equation  (20)  gives,  for 
p  ■  0, 

o  =  ±  33.1° 


which  is  the  interblade  phase  angle  for  resonance  found  by  Verdon 
and  Caspar. 

Results  for  the  flat  plate  case,  computed  using  Equation 

(24)  and  100  blades,  are  shown  in  Figure  3.  In  this  calculation 

the  total  number  of  blades  in  the  numerical  computation  is  13  and  m 

3.  A  transonic  case  for  a  NACA0012  section,  M  =  0.75,  h  =  2.0 

00 

is  shown  in  Figure  4. 

Nonunique  Solutions 


Over  the  last  six  years  nonunique  solutions  to  the  transonic 
potential  equation  have  been  found.  (See  References  11  and 
12).  In  the  initial  cases  these  consisted  of  lifting  solutions 
for  the  flow  around  symmetric  airfoils  at  zero  angle  of  attack. 
Later,  nonunique  solutions  were  found  for  nominally  lifting 
conditions.  These  nonunique  solutions  have  been  found  for 
isolated  airfoils  and  for  airfoils  in  wind  tunnels;  unsteady 
solutions  have  also  been  reported  in  Reference  13.  In  the 
present  work,  nonunique  solutions  have  been  found  for  cascade 
solutions. 

In  order  to  find  the  nonunique  solutions  the  cascade  code, 
CTRAN,  developed  under  previous  work  on  the  contract,  was  used. 
Also,  a  variant  of  the  code  XTRAN2L,  modified  to  give  periodic 
boundary  conditions,  was  used.  Only  steady  flows  are  considered 
at  the  moment.  For  a  gap  of  50  chord  lengths  the  range  of  Mach 
numbers  for  a  nonunique  solution  are  shown  in  Figure  5.  In 
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; 


15).  The  lift  induced  on  a  blade  is  proportional  to  the  induced 
upwash  parameter,  whj^h  is  given  by  Equation  (16).  Hence  it 
follows  that 


CT  =  cr  ISLI  l/2eiK(  |n|-m)h  (17 

L  L  i  n 1 

n  m 

where  the  relation  y  =  nH  is  used,  and  m  is  a  given  blade  beyond 

which  Equation  (16)  is  .a  good  approximation.  It  is  assumed  that 

CT  is  known.  Note  that  the  actual  value  of  v(x,y)  is  not  given 
Hn 

by  Equation  (16)  since  this  result  only  holds  if  there  are  no 
other  blades.  However,  the  actual  upwash  will  vary  with  y  in  a 
form  similar  to  Equation  (16).  Using  Equations  (14),  (15)  and 
(24)  gives 


m-l  '  00 

CT  =  l  CT  e-ina  +  l  CT  (£)1//2  exp  [  iK ( n-m)K-ino] 
n=-m+l  n  n=m  m 

+  l  Cr  (S) 1/2exp[-iK(n+m)h-ino)  (18) 

n=-«  -m  n 

It  can  be  seen  that  the  effect  of  the  cascade  blades 
diminishes  as  n-^/2  and  hence  a  large  number  of  blades  is 
necessary  in  the  initial  computation  of  CL  unless  the 
approximation  of  Equation  (17)  is  used.  It  is  also  of  interest 
to  examine  further  Equation  (1C).  The  two  infinite  sums  are 
bounded  (see  Appendix  4)  unless 

±  kE  -  o  =  2pn  (19) 

where  p  is  an  integer.  Thus  for  an  interblade  phase  angle,  given 
by 


o  *  ±  kK  -  2pn 


(20) 
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cascades  are  described.  The  results  did  not  agree  very  well  with 
those  of  Verdon  and  Caspar  (Reference  9)  and  hence,  some 
examination  into  the  causes  of  the  discrepancy  is  necessary. 


In  Reference  10  the  computed  data  are  for  a  seven  blade 
cascade  and  it  is  possible  that  this  restricted  number  of  blades 
is  a  cause  of  the  discrepancy.  This  aspect  is  elaborated  below. 

The  lift  given  by  the  theory  in  Reference  1  is 
00 

C  (t)  *  l  CL  (t  -  n  a)  (14: 

n=-«  n 

where  CT  is  the  lift  induced  at  blade  zero  by  the  motion  of  the 
Ln 

ntfl  blade  with  the  other  blades  stationary  and  o  the  interblade 
phase  angle. 


If  the  flow  is  harmonic  in  time. 


CL(t)  =  CLeJ 


CT  (t  -  no)  =  CT  e 
A-»_ 

n  n 


iw(t  -  no) 


In  a  cascade  the  unsteady  flow  induced  on  blade  zero  by  the  ntn 
blade  can  be  characterized  by  the  asymptotic  form  of  the  unsteady 
flow  provided  the  nfc^  blade  is  sufficiently  far  from  blade  zero 
(for  example  several  chord  lengths).  The  asymptotic  behavior  for 
the  induced  normal  velocity  is  (Appendix  3) 


v(x,y)  ~  A  — 1/2  , 
(K|y|1/2) 

—  2  1/2 

where  A  is  a  complex  constant  and  y  =  (1-M  )  y 


iK  y 


This  result 


is  obtained  from  classic  subsonic  airfoil  theory  (see  Reference 


Steady  pressure  distribution  through  a 
cascade  =  0.845,  a  =  0,  h  =  30. 

Figure  1 
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Unsteady  pressure  distribution  through  a 

cascade  M  =  0.845,  a,  =  0.25°,  h  =  30. 

00  1 

Figure  2a 


Variation  of  nonunique  behavior  with  Mach 


Variation  of  nonunique  behavior  with  gap 


0.6 


30 


Unsteady  nonunique  behavior 
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APPENDIX 


1 .  INTRODUCTION 


The  most  common  methods  of  predicting  steady  flow  aero¬ 
dynamic  characteristics  at  transonic  speeds  are  either  the 
Transonic  Small  Disturbance  (TSD)  theory  (Ref.  1)  or  the 
Full  Potential  Equation  (FPE)  theory  (Ref.  2) .  The  more 
accurate  Euler  equations  solutions  (Ref.  3)  are  expensive  to 
obtain,  although  for  flows  with  strong  shock  waves  such  solu¬ 
tions  axe  essential.  The  FPE  theory  is  based  on  the  assumption 
that  the  flow  is  isentropic  and  irrotational  and  generally  has 
a  (numerically)  exact  treatment  of  wing  boundary  conditions. 

The  TSD  theory  is  an  approximation  to  the  FPE  theory  and  thin 
wing  boundary  conditions  axe  used  in  the  solution  procedure. 

One  of  the  advantages  of  the  TSD  theory  is  the  flexibility  in 
deriving  the  approximate  equation.  This  flexibility  is  generally 
utilized  by  a  choice  of  a  transonic  scale  parameter. 

The  basic  assumptions  of  isentropy  and  irrotationality  in 
both  these  theories  are  only  valid  when  there  are  no  shock 
waves  (Ref.  4)  in  the  flow.  Although  the  potential  theory  is 
valid  strictly  only  for  shock-free  flows,  the  results  obtained 
by  such  a  theory  for  flows  with  shock  waves  are  sufficiently 
close  to  experimental  data  for  practical  use,  provided  the  shock 
waves  are  weak.  The  generally  accepted  definition  of  a  weak 
shock  is  when  the  local  Mach  number  just  ahead  of  the  shock  is 
less  than  1.3.  Thus,  when  both  TSD  and  FPE  solutions  are 
compared  to  the  more  realistic  Euler  equation  solutions  it  is 
found  that  the  agreement  is  satisfactory  provided  that  the 
basic  restriction  to  weak  shock  waves  is  not  violated.  The 
use  of  thin  wing  boundary  conditions  can  also  introduce  errors 
into  the  TSD  solutions. 


If  the  flow  has  strong  shock  waves,  however,  then  there  is 
considerable  disagreement  between  both  TSD  and  FPE  solutions 
and  Euler  equations  solutions.  Generally  the  predicted  shock 
location  for  the  potential  theories  is  much  further  aft  than 
that  for  the  Euler  equation  solutions.  The  causes  of  the 
error  in  the  shock  location  in  the  steady  TSD  theory  for  two- 
dimensional  flow  have  been  examined  in  Reference  4  where  the 
basic  justification  for  a  correction  procedure  has  been  derived. 
The  basic  hypothesis  of  the  theory  is  that  the  error  in  shock 
location  is  primarily  due  to  the  stronger  shock  predicted  by 
TSD  theory  compared  to  the  shock  strength  of  the  Euler  equations. 
It  is  also  assumed  that  if  the  shock  strength  is  suitably 
corrected  then  the  shock  location  should  be  approximately 
correct.  These  correction  theories  have  been  applied  to  both 
steady  (Ref.  5)  and  unsteady  (Ref.  6)  solutions  of  the  TSD 
equation  and  to  the  steady  FPE. 

The  commonly  used  potential  theories  are  based  on  Crocco's 
theorem  which  states  that  if  the  entropy  gradient  in  the  flow 
is  negligible  then  an  inviscid  flow  exists  with  negligible 
vorticity  and  hence  a  velocity  potential  derived  under  the 
assumption  of  zero  vorticity  gives  a  good  approximation  to 
the  flow.  However,  application  of  Crocco's  therem  requires  that 
mass,  momentum  and  energy  be  conserved  and  since  in  a  potential 
solution  with  shock  waves  momentum  is  not  conserved,  it  can 
be  seen  that  the  transonic  potential  theory  is  not  consistent 
with  its  basic  assumptions.  However,  since  transonic  potential 
theory  does  give  good  results  for  flows  with  weak  shock  waves 
its  use  can  be  justified  on  the  basis  of  practicability  alone. 

A  more  flexible  criteria  for  developing  consistent  theories  is 
to  introduce  (Ref.  4)  the  idea  of  minimizing  a  weighted  com¬ 
bination  of  source  errors  at  the  shock  wave.  For  example,  if 
mass  and  energy  are  conserved  through  a  shock  then  the  conven¬ 
tional  potential  theory  results,  or,  if  a  suitable  weighted  sum 
of  the  mass,  momentum,  and  energy  errors  are  put  to  zero  then  a 
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potential  equation  with  a  specified  shock  jump  results.  Thus, 
a  potential  theory  with  a  Rankine-Hugoniot  shock  jump  can  be 
obtained . 

In  Reference  5  the  correction  to  the  TSD  equation  is 
obtained  by  computing  two  steady  state  solutions  and  then 
using  an  interpolation  technique  to  give  the  required 
solution.  This  technique  is  not  really  feasible  for  unsteady 
flow  since  the  correction  procedure  is  required  for  each 
time  step  in  the  two  TSD  solutions,  with  different  scaling 
parameters,  and  an  interpolation  scheme  derived  for  discon- 
tinous  transonic  flows.  Examples  of  steady  flows  with 
strong  shocks  computed  with  this  method  agree  satisfactorily 
with  the  Euler  equation  solutions,  although  the  use  of  the 
thin  airfoil  boundary  conditions  in  the  TSD  theory  can  give  rise 
to  errors  near  the  leading  edge . 

An  extension  of  the  basic  correction  procedure  for  unsteady 
transonic  flow  is  given  in  Reference  6  where  the  low  frequency 
theory  of  Ballhaus  and  Gocrjian'  {Ref.  8)  is  extended  by  adding 
a  formally  negligible  third  order  term  to  the  theory.  The 
improvement  in  the  potential  theory  for  flows  with  strong 
shock  waves  is  remarkable.  The  conventional  theory  failed  to 
give  an  answer.  However,  this  theory  still  has  some  problems, 
in  particular,  the  agreement  of  the  pressure  distributions  with 
the  results  of  the  Euler  equations  is  unsatisfactory  ahead  of 
the  shock  wave.  This  is  attributed  to  the  global  nature  of  the 
strong  shock  correction  affecting  the  flow  ahead  of  the  shock 
wave  as  well  as  weakening  the  conventional  potential  shock. 

Although  all  of  the  above  theories  have  been  vindicated 
by  computations  the  number  of  examples  computed  is  very  small. 
Before  any  of  these  theories  can  be  used  in  practice  with 
some  confidence  a  much  more  detailed  validation  is  required. 


In  the  present  work,  therefore,  the  application  of  the  strong 
shock  theory  to  unsteady  flows  will  be  tested  in  much  greater 
detail  than  has  been  done  at  present. 

In  the  present  study  it  is  found  that  the  pressure  error 
in  the  previous  study  is  due  to  one  of  the  additional  parameters 
being  unconstrained.  Also,  it  has  been  found  that  the  theory 
of  Reference  7  is  not  easily  extended  to  unsteady  flows  (see  the  Appendix)  . 
The  theory  is  extended  to  treat  high  frequency  flows.  The  results 
of  the  modified  theory  show  a  marked  improvement  over  the 
earlier  work  of  Reference  6. 


2.  BASIC  EQUATIONS 

The  basic  differential  equation  used  in  the  present 
analysis  is 


(a  +  b<fc  +  c<p  )  d  +  d  +  Ad  4.  +  Bd,.*  =  0 
1  yx  x  xx  Yyy  rxt  tt 


where  a,  b,  c  are  parameters  to  be  chosen, 


A  =  -2Kfv/c2/3 


B  =  -M2v2/62/3 

CO  ’ 


where  M  is  the  free-stream  Mach  number  and  v  is  the  reduced 

cc 

frequency.  If  6  is  the  airfoil  thickness  parameter,  c  is  the 
airfoil  chord  and  Uro  the  free-stream  velocity,  then  x,y,t,d  are 
normalized  with  respect  to  c,  c/5^3,  w  3 ,  c52/,3Ud,  where  u 
is  the  oscillation  frequency. 


The  high  frequency  boundary  condition  is 

c  (x,iO,t)  =  P  (x,±0,t)+^r-  (x,±0,t) 
y  t  x  ct 


(3) 


where  v  =  f(x,tO,t)  denotes  the  upper  and  lower  surfaces  of  the 
airfoil  in  motion,  respectively.  The  far  field  boundary  condi¬ 
tion  is  a  nonref lective  boundary  condition  of  the  type  developed 
by  Kwak  (Ref.  9)  and  is  further  developed  by  Whitlow  (Ref.  10). 
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In  Figure  (2) ,  the  steady  pressure  distribution  and  the 
unsteady  pressure  jump,  <Cp-C")  -  (C*  -  C~)  steady  is  shown  for 
a  NACA  0012  airfoil  at  Mw  =  0.8,  ac  =  0.7°  with  =  0.25°  and 
v  -  0.6.  It  may  be  seen  that  the  steady  state  shock  has  moved 
location  considerably.  There  is  a  considerable  change  in  the 
unsteady  pressure  jump. 

Figure  (3)  shows  the  steady  and  unsteady  results  for  a 
NACA  64  A006  at  =  0.875  at  dQ  =  0.0  with  a  flap  deflection 
of  1.0°  and  v  =  0.470;  the  flap  hinge  is  at  25%  of  chord. 

This  is  case  10  from  Ref.  11.  It  can  be  seen  that  while  there 
is  a  relatively  small  change  in  the  steady  pressure  distri¬ 
bution,  there  is  a  considerable  effect  in  the  unsteady  result, 
mainly  due  to  the  change  in  location  of  the  shock. 

In  Figure  (4) ,  the  results  are  shown  for  the  flow  around 
a  NACA  64A010  airfoil  at  =  0.796,  =  0.0  with  a  pitching 

amplitude  of  1.01°  and  v  =  0.404.  The  results  are  similar  to 
those  of  the  previous  example.  This  is  case  6  of  Reference  11. 

Finally,  in  Figures  (5)  and  (6) ,  the  results  for  an 
NLP.  7301  airfoil  at  M„.  =  0.721  and  ac  =  -0.19°  are  shown. 

These  are  cases  8  and  13  frcm  Reference  11.  There  is  a 
double  shock  in  this  example,  the  first  of  which  moves  for¬ 
ward  due  to  the  present  modifications  of  the  method  for  steady 
flow,  although  the  unsteady  results  do  not  differ  significantly. 
In  Figure  (5),  the  airfoil  is  oscillating  in  pitch  at  v  =  0.362 
with  amplitude  0.5°.  The  unsteady  result  differs  considerably 
from  the  original  result.  In  Figure  (6) ,  there  is  a  flap  hinged 
at  75%  of  chord  oscillating  at  v  =  0.362  with  amplitude  1.0". 

It  should  be  noted  that  this  result  is  for  the  third  cycle;  in 
the  modified  method  the  fourth  cycle  civerged--the  reasons  are 
not  known. 
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4.2  Unsteady  Flow 

The  code  uses  the  theory  outlined  elsewhere  except 
that  the  values  of  a  and  8  are  computed  only  on  the  upper  and 
lower  surfaces  of  the  airfoil  and  these  are  used  throughout 
the  appropriate  half  plane.  This  device  is  for  coding  simplicity 
and  assumes  that  the  unsteady  effects  on  the  shock  wave  can  be 
characterized  by  the  behavior  at  the  shock  foot.  It  is  also 
possible  that  instabilities  in  the  algorithm  could  arise  if 
a, 8  were  varied  along  the  shock.  An  instability  did  arise  in 
the  steady  flow  part  due  to  varying  b,c  in  an  injudicious 
manner . 

5.  DISCUSSION  OF  RESULTS 

One  of  the  main  purposes  of  the  present  work  is  to  test 
the  ideas  of  the  present  theory  over  a  range  of  cases.  These 
cases  include  some  from  Reference  (10)  and  compare  the  present 
method  with  the  unmodified  method. 

The  first  case  is  the  steady  flow  around  a  NACA  0012 
section  at  =  0.8,  a  =  1.25°,  a  case  for  which  Euler  equation 
solutions  are  available.  The  comparison  is  shown  in 
Figure  (1) .  It  can  be  seen  that  the  present  result  improves  the 
agreement  with  the  Euler  solution  by  weakening  the  upper  surface 
shock  strength  and  moving  its  location  forward  slightly  ahead  of 
the  Euler  solution.  A  similar  overcorrection  is  seen  on  the 
lower  surface.  This  may  be  due  to  the  strength  being  below 
the  "cut  off"  strength,  in  this  case  |  |  =  5,  and  hence 

the  present  correction  is  not  applied  to  the  lower  surface  flow. 

It  should  also  be  noted  that  this  solution  is  not  fully  converged, 
the  residual  being  about  S  x  10  however,  the  solution  is  not 
obviously  diverging.  It  should  also  be  noted  that  the  pressure 
overshoot  shown  in  the  earlier  work,  Ref.  6,  has  been  similarly 
overcorrected. 


4.  COMPUTATIONAL  PROCEDURE 

The  theory  developed  in  the  preceding  section  was  implemented  in 
the  computer  code  XTRAN2L  (Ref.  10),  as  follows: 


4.1  Steady  Flow 

When  a  shock  appears  the  strength  of  the  small  disturbance 
shock  is  computed.  If  |  $3  ~  $31  is  less  than  a  specified  value 
then  no  alteration  is  required;  at  present  the  specified  value 
is  0.1.  If  the  shock  is  sufficiently  large,  then  the  constant  c 
is  chosen  to  give  a  shock  with  the  Rankine-Hugoniot  shock 
strength.  This  is  accomplished  using  the  formulae  of 
Equation  (16) .  Incidently,  in  some  of  the  examples  the  exact 
nonlinear  formula  for  Mg  gives  an  imaginary  value;  the  code  was 
modified  to  use  the  linear  equivalent  formula 

.,2  ,  .  .2/3, 

=  M  +  kc  7  $ 
e  00  x 

where  k  is  the  transonic  parameter.  If  c  is  found  to  exceed 
, b;/6,  then  c  is  kept  at  its  previous  value  and  b  is  modified 
to  achieve  the  correct  shock  strength.  This  reduces  the 
likelihood  of  multiple  parabolic  points. 


In  order  to  retain  a  stable  algorithm  the  coefficients 
b,c  are  updated  only  every  50  iterations.  Furthermore,  the 
values  on  the  upper  and  lower  surfaces  are  used  throughout  the 
respective  half  plane  until  the  residual  is  five  times  the 
criteria  for  convergence  or  that  the  number  of  iterations 
exceeds  90%  of  the  total  maximum  specified  iterations.  The 
correct  field  values  are  then  used  to  increase  accuracy  and  to 
allow  the  intermediate  field  solution  to  reduce  to  its  classic 
form,  so  that  the  boundary  condition  treatment  is  not  compromised. 

For  each  sweep  over  the  airfoil  the  coefficients  are  updated 
only  if  the  shock  strength  exceeds  by  10%  the  strength  of  a 
previous  shock  calculated  on  the  same  sweep.  This  allows  for 
multiple  shock  waves  and  updates  the  coefficients  only  for  the 
first  shock  or  the  subsequent  shocks  if  they  are  stronger  by  10% 


he  first. 
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f (u)  =  a  +  bu  +  cu^ 

then 

3  f 

r-r  <  0  for  a  range  of  u 

d  u 

or 

b  +  2cu  <  0 

Let  umax  be  the  maximum  value  of  u  in  the  range,  then 
-b  .  . 

c  2u -  C  Posltlve  (30) 

max 

if  umax  =  4'  say'  then 

jff  i  !/e  (3D 


In  practice  a,b  are  approximately  the  same  order  of  magnitude 
and  hence  the  inequality  of  Equation  (29)  is  the  more  restric¬ 
tive.  Thus, 

2 

c  <<  b  /a 


or 


max 


1 

8 


a 


(32) 


which  is  used  in  the  analysis;  the  value  "1/8"  is  arbitrary. 
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3.  CONSTRAINTS  ON  THE  COEFFICIENTS 


There  is  a  danger  of  the  quadratic  in  Equation  (1)  term 
having  multiple  roots  that  are  physically  not  unreasonable  but 
which  could  destroy  the  algorithm  stability. 

The  parabolic  points  of  Equation  (1)  are 


u  =  _0i 


-  4ac] 
2c 


<<  1  and  a  -  £(1)  then  a  Taylor's  expansion  gives 


«•-! 


c  b 


The  proper  solution  is  the  first  root  since  it  has  the  correct 
behavior  as  c  +  0. 

In  order  to  avoid  the  second  root  it  must  be  much  larger 
or  smaller  than  the  first  so  that  it  is  not  likely  to  appear  in 
a  physically  reasonable  situation.  In  general,  a  >  0,  b  <  0 
and  | a]  and  jb|  are  of  similar  magnitude  and  thus 

if  -  -  >  0  then  -  -  >>  -  £ 

C  CD 

and 

if  -  -  <  0  then  -  -  <<  0 

c  c 


Hence  for  the  second  root  to  be  avoided,  the  following 
constraints  must  be  met. 

jcl  «  |bl  c  negative 

2  (29) 

c  «  —  c  positive 

A  further  constraint  is  that  the  nonlinear  term  should  be 
monotonic  through  a  range  of  u.  Thus,  if 


If  these  equations  are  satisfied  then,  to  first  order 
in  shock  motion,  the  shock  will  have  a  strength  equal  to  the 
Rankine-Hugoniot  value. 

For  steady  flow  a  is  kept  equal  to  its  traditional  value, 
and  b  and  c  are  altered  to  satisfy  equation  (16)  .  There  are 
some  constraints  on  b,c  and  these  are  discussed  later. 

For  unsteady  flow,  c  is  kept  equal  to  its  steady  state 
value  and 


b  =  bs(*W  >°  +  exs 


where  the  subscript  s  denotes  a  steady  state  value.  The  a  and 
g  are  found  by  satisfying  Equations  (20)  and  (21) .  Thus, 

+  f  r  i  f3aEi  i  r  + 

aVV  |bs  1  "  2  77+  +  C  2<xs  •  CES 

s  36  s 

l  sJ 


— 


K. 


/U  -  0Es/2) 


[bs/2  + 


o  r  + 

4  C0„  1  -r*  +  A  /u*  ~  0  /2) 

3  Es  3xs  xs  Es 


where 


ScE  es 

-  =  -vc_  —5 -  +  1 

ox.  ES  M2  -  1 


where  M  is  the  steady  state  Mach  number  ahead  of  the  shock. 
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(Y+1)M„ 


0  =  ~  ^  =  17^7?  (M*  •  2T  +  —  cpJ}/{(1  +  vV^> 

(15) 

The  problem  is  now  reduced  to  replacing  a  in  Equation  (8)  with 


Following  the  philosphy  of  the  earlier  work,  the  steady 
shock  strength  is  fixed  first.  Thus, 

a  +  b(4>+  -  o£s/2)  +  c(<j>*  -  ^  o£s  +  °£S/3)  «  0  ' 

s  S  S 

This  allows  the  constants,  a,  b,  c  to  be  chosen. 

Next,  the  shock  strength  is  expanded  about  its  steady 
state  value.  Thus,  if 


F  Ux,o,xs)  =  0 


denotes  Equation  (8) ,  then 
F(^,c,xs)  =  F(<^  ,  Oj 


Es' 

0)  +  xs 

[cxs 

<}>+  ) 
V 

?F  + 

~r  3c_- 

cF  E 

^ a. 

cF  c. 


E  oXS 


or,  since 


F(<  '  °Es '  0)  =  ° 


?F  +  _cF_ 

ScE  SXs 


5  F  +  _c_F_  ffE 
3<  5°E  ^ 
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where  Me  is  Mach  number  just  upstream  of  the  shock  and  is 
defined  as 

M  =(U+  -  x  )/a+  (10) 

c  S 

•  +4. 

where  x  is  the  physical  shock  speed,  U  and  a  are  the  velocity 
s 

and  speed  of  sound  just  ahead  of  the  shock,  respectively.  Now 
mJ  =  (U+  -  xg)2/{  af[l  +  lli  m2(1  -  U+2/vl )]}  (11) 

or 

M2  =  M2  (1  -  x _/U+)2  (12) 

e  es  s 

where  M„  is  the  steady  Mach  number  with  the  steady  value  of  U+ 
es 

replaced  by  its  unsteady  value. 

In  the  scaled  coordinate  of  the  present  problem,  the 
physical  shock  speed  is  replaced  by 

viu,,, 

S  00 

Hence , 

M2  =  M2  [1  -  vx  / (U+/U  ) ) 2  (13) 

e  e  _  s  00 

s 

Returning  now  to  Equation  (9);  if  the  small  disturbance  pressure 
relation, 

C  =  -2  ($*  +  2v<t>. )  o2/3  (14) 

P  X  w 

is  used,  then,  with  the  help  of  Equation  (8) 
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The  shock  jump  relation  is 


U[*J  +  §1**3  +  |  “  [*yl 


ax 


{A($x]  +  B  [  1 3  } 


dx 

at 


=  0 


(4) 


where  [  ]  denotes  a  jump  across  the  shock.  Since  (dx/dt)  =  x  , 

5  S 

the  shock  speed,  is  continuous  through  the  shock,  then 


(5) 


and  hence  for  a  normal  shock  Equation  (4)  becomes 

{£  +  I  U*  +  $">  +  f  +  *x6x  +  ^x2)}  '  *s{A  +  BV 


0 

(6) 


If  a  shock  strength,  c,  is  defined  by 


c  =  e  -  t 
x  x 


(7) 


Tier.  location  (6)  becomes 


{a  +  b(£  +  -  c/2)  +  c  {$+  -  +  o2/3)) 

X  .  X  A 

-  x  {A  +  Bx  )  =  0 
s  s 


:8) 


Consider  now  the  one-dimensional  Euler  shock  strength.  The 
pressure  jump  is  given  by 


1) 


(9) 
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APPENDIX 


NOTES  ON  THE  EXTENSION  OF  THE  KLOPFER-NIXON 
METHOD  TO  UNSTEADY  FLOW 


In  Reference  7  a  method  is  developed  to  modify  the  steady 
full  potential  equation  to  treat  strong  shock  waves.  The 
results  of  this  method  agreed  better  with  solutions  of  the 
Euler  equations  than  the  steady  version  of  Reference  (6)  and 
it  is  therefore  instructive  to  see  whether  some  of  these 
ideas  can  be  extended  to  time  dependent  flows. 

Consider  the  two  dimensional  Euler  equations: 

Sp  .  _  0 

St  Sx^ 

Spu.  3pu.u.  _  3o 

St  cXj  Sx± 

il  +  u.  iH_  ,  o 

St  1  oxi 

v:here 


H  =  h  +  j  uiui 


Y  -  1 


For  steady  isentropic  flow  the  energy  equation  can  be  integrated 
along  a  streamline  to  give 


h  +  2  uiui  =  constant  =  hQ 


and  the  isentropic  relation 
p/pY  =  constant  =  k  p^/p'1 


car.  be  used  to  give 


,  |  y  -  1  „2  ,  UiUi 

c  p  1  +  j — Ml-  j— 


2 - »  *  2  I 

q=c  j, 


Y-l  1 
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This  relation  is  then  used  in  conjunction  with  the  conserva¬ 
tion  of  mass  equation  and  the  irrotationality  equation 

ui  =  7x7  '  <*> 

To  give  an  equation  for  $ .  In  the  usual  case  k  =  1  but  for  the 
strong  shock  case,  k  is  allowed  to  jump  at  the  shock  wave. 

For  unsteady  flow  the  energy  equation  cannot  be  integrated 
and  the  necessary  density/velocity  relation  is  obtained  from 
the  momentum  equation. 

Consider  the  x  momentum  equation  in  two  dimensions 


3u  +  u3u  v3u  :  -  I  IE 
St  3x  3y  p  3x 


Using  irrotationality  and  the  isentropic  relation  gives 


4r5(“2  +  v2)+rn  [k(p/pj^  o 


If  k  is  constant,  then  Equation  (10)  reduces  to  the  form 


33612  2  o  I 

It T  +  i  (u  +  v  )  +  - —  k  -f- 

cx  (y-1)p.  lP“j 


which  can  be  integrated  to  give  the  unsteady  Bernoulli  equation. 
If  k  is  allowed  to  jump  through  the  shock  then  the  equation  is 


—  ,  fL£  +  1  (u2  +  v2)  +  — —  k 

ix  3 1  2  '  Y-l  P~ 


Y-J-  P»  P« 


+  ^  -B-  =  0 

p  p  oX 

cc  Kco 
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which  does  not  give  a  simple  density/velocity  relation.  Thus, 
the  application  of  the  theory  of  Ref.  7  to  unsteady  flow  is 
more  difficult  than  first  envisaged. 
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In  recent  years  multiple  solutions  to  the  numerical 
approximation  to  the  full  potential  equations  have  appeared  in 
the  literature  (Refs.  1  and  2).  Initially  the  phenomena  appeared 
in  computations  of  the  flow  over  a  symmetric  airfoil  at  zero 
angles  of  attack  whre  two  lifting  solutions  were  present  in 
addition  to  the  expected  nonlifting  solution.  In  Reference  2 
some  results,  for  a  nonsymmetric  airfoil,  a  RAE  2B22  section,  are 
also  presented.  Steinhoff  and  Jameson  (Ref.  1)  suggested  that 
the  change  from  one  of  the  solutions  to  another  is  discontinuous 
and  noted  a  hysteresis  effect  indicating  that  the  lift 
coefficient  (CL)  depended  on  whether  the  angle  of  attack  (a)  was 
increasing  or  decreasing.  More  recent  work  is  by  Salas  (Ref.  2) 
who  has  extended  the  computations  of  the  flows  considered  by 
Steinhoff  and  Jameson  (Refs.  1)  to  show  that  it  is  possible  to 
construct  a  smooth  CL  -  a  curve  connecting  the  three  solutions 
for  a  symmetric  airfoil. 

The  investigations  noted  above  are  meticulously  performed 
and  are  essentially  numerical  experiments.  There  is  a  limited 
amount  of  understanding  that  can  be  gained  from  such  experiments 
and  consequently  a  more  analytic  technique  may  yield  additional 
information.  Furthermore,  although  the  numerical  results  are 
invaluable  they  do  not  exclude  the  possibility  that  the  multiple 
solutions  are  due  to  the  numerical  approximation  to  the 
differential  equation.  The  present  investigation  is  based  on  the 
integral  equation  formulation  (Ref.  3)  which  allows  some  degree 
of  insight  into  the  problem. 

The  objective  of  this  paper  is  to  suggest  a  possible  reason 
for  the  multiple  solutions  based  on  the  existence  of  an 
eiqensolut ion  in  the  transonic  small  disturbance  ( TSD ) 
equation.  It  is  shown  that  a  ficticious  lift  can  be  added  to  a 
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"real"  solution  to  give  a  multiple  of  solutions.  It  is  also 
indicated  that  the  methods  used  to  stabilize  a  numerical  solution 
are  likely  to  indicate  the  appearance  of  the  eigensolut ion . 

2.  INTEGRAL  EOUATION  ANALYSIS 

The  transonic  integral  eguation  method  of  Reference  3  is 
only  applicable  to  the  transonic  small  disturbance  (TSD)  eguation 
rather  than  the  full  potential  equation  (FPE)  that  is  used  in  the 
earlier  work.  Consequently  the  first  step  is  to  reproduce 
multiple  solutions  using  the  TSD  equation.  Once  these  solutions 
are  obtained  they  can  be  analyzed  using  the  ideas  of  the 
transonic  integral  equation  theory. 

2.1  Multiple  Solutions  for  Small  Disturbance  Theory 

Since  it  is  easiest  at  present  to  use  the  inteqral  equation 
theory  to  analyze  small  disturbance  theory  it  is  necessary  to 
reproduce  the  multiple  solutions  using  the  TSD  equation.  This  is 
achieved  by  usina  the  computer  code  TSFOIL  (Ref.  4)  which  solves 
the  TSD  equation  using  the  conservative  Murman-Cole  algorithm  and 
grid  sequencing.  The  multiple  solutions  are  found  for  a 
symmetric  airfoil  at  zero  angle  of  attack  by  imposing  a  1°  angle 
of  attack  on  the  coarse  grid  solution  and  then  putting  the  angle 
of  attack  equal  to  zero  in  the  medium  and  fine  arid  operations. 

It  is  found  that  such  a  device  leads  to  multiple  solutions  over  a 
small  ranae  of  Mach  numbers.  Such  solutions  have  been  found  for 
a  11.8%  Joukowski  airfoil  and  a  NACA  0012  airfoil.  As  a  test  for 
convergence,  the  solution  for  the  Joukowski  airfoil  at 

M  =  0.85  was  converged  to  a  residual  of  10”^.  Krupp  scaling 
00 

was  used  in  these  results;  the  default  Krupp  grid  is  used,  which 
has  grid  dimensions  of  77  x  56.  An  example  of  a  multiple 
solution  is  given  in  Figure  1. 
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2.2  Transonic  Inteqral  Equation  Theory 


Since  the  analysis  of  the  multiple  solutions  is  based  on  the 
integral  equation  theory  it  is  helpful  to  outline  the 
f ormu la t ion . 


For  Krupp  scaling  the  TSD  equation  can  be  written  as 


( 1  —  <(>  )  $  +  =0 
x  xx  yy 


(1) 


where  if  $  is  the  perturbation  velocity  potential  in  a  Cartesian 
coordinate  system  (x,y)  then 


where 


<t>  =  k/B 
x  =  x 
y  =  3v 


2  - 


(2) 


k  = 
v  = 


( ,  + 

(  ?3  7k  )  v 


1  .  75 


(3) 


If  U  =  y 

X 

components 


v  =  ;  then  the  physical  perturbation  velocity 


u,v)  are  given  by 


u  =  ( 3  /k  )u 
V  =  ( e3/k  )v 

In  the  formulation  of  Equation  (1)  the  sonic  line  is  given  by 


u  =  4>  =  1 
x 


The  boundary  conditions  for  Equation  (1)  are  that 

?  ,  <p  *  0  on  the  far  field,  that  the  tanqency  condition 
x  y 


(4) 


(5) 


$  ( x , ±0)  =  y^ ( x , ±0)  -  A 


(6) 
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It  is  of  interest  to  note  that  if  a  stabilizina  tern  of  the 
type  used  in  TSFOIL,  namely  £4>xt  is  added  to  Equation  (1),  where 
t  is  an  artificial  iteration  time  and  e  is  a  parameter,  then  an 
analysis  similar  to  that  given  in  this  section  shows  that  this 
term  will  assist  the  formation  of  a  nonunique  solution.  This  is 
possibly  the  reason  why  the  real  solution  cannot  sometimes  be 
computed.  The  same  effect  is  also  present  if  the  dissipation  due 
to  the  truncation  error  in  the  upwind  difference  scheme  is 
indicated  in  the  analysis.  In  other  words,  a  stable  conservative 
algorithm  is  likely  to  initiate  the  appearance  of  the 
eigensolution  g  in  the  solution. 

3.  SUMMARY 

The  TSD  equation  can  admit  eigensolutions  that  satisfy  all 
of  the  boundary  conditions  generally  found  in  such  problems; 
there  are  an  infinite  number  of  these  eigensolutions.  For  a 
nonunigue  solution  to  exist  certain  consistency  re  j  i  remen  ts  must 
be  met.  These  eigensolutions  provide  lift  and  act  like  an 
additional  asymmetric  source  term  in  the  flow  field.  If  real 
lift  is  present  the  fictitious  component  provided  by  the 
nonuniqueness  appears  as  a  simple  additive  term  to  the  real 
lifting  component;  there  is  no  way  to  distinguish  or  uncouple  the 
two  components  in  a  given  numerical  solution. 

If,  in  a  supercritical  flow,  an  eiqensolution  does  appear, 
the  location  of  the  shock  waves  changes  from  their  real 
location.  The  nonuniqueness  may  be  removed  in  some  steady  cases 
by  usinq  a  nonconservative/conservative  algorithm  which  may 
inhibit  the  appearance  of  such  solutions.  However,  there  is  no 
guarantee  that  this  will  work  in  all  cases.  This  partial  cure 
only  works  for  steady  flows;  it  is  possible  that  for  unsteady 
flows  that  a  "strong  shock"  theory  of  the  type  advocated  in 
Reference  8  could  remove  any  nonuniqueness.  The  (necessary) 
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lower  surface  shock  moves  forward.  It  follows  that  the  source 
term  introduced  by  the  shock  on  the  lower  surface  now  acts  at  a 
further  forward  location  while  the  countering  source  on  the  upper 
surface  moves  aft.  Consequently/  there  is  a  considerable  region 
between  the  shocks  for  which 

A [ o . H( x-x  ) ]  <0 

i  Si 

This  counters  the  positive  g  and  reduces  the  tendency  of  the 
nonunigue  solution  to  appear.  However,  it  is  necessary  to  point 
out  that  the  nonconservative  source  error  must  be  large  enough  to 
counter  the  g  terms.  The  nonconservative  algorithm  only  inhibits 
a  nonunique  solution  from  appearing;  it  will  not  necessarily 
remove  an  existing  error.  A  similar  analysis  can  be  performed  if 
the  fictitious  lift  is  negative.  This  hypothesis  was  tested  by 
using  the  computer  code  TSFOIL  (Ref.  4)  which  has  both 
conservative  and  nonconservative  algorithms.  A  composite 
alaorithm,  given  by 

Algorithm  =  X  (conservative)  +  (1  -  X)  (nonconservative)  (65) 

was  used.  The  parameter,  X,  was  taken  to  be  given  by 

X  =  -e|[CL(t+At)  -  CL(t)]|  +  1  (66) 

where  t  is  the  artificial  iteration  time  and  At  is  the  iteration 
step.  At  convergence  X  =  1  and  the  solution  is  conservative. 

It  was  found  that  for  the  nonuniqut  solutions  to  be  avoided  in  a 
computation  of  the  example  in  Figure  1  efC^tt+tAt)  =  (^(t)]  had 
to  be  of  order  unity  during  most  of  iteration.  Smaller  values 
did  not  inhibit  nonuniqueness.  Hence  it  is  suggested  that  the 
"classic"  deqree  of  nonconservative  algorithm  is  necessary  to 
stop  the  appearance  of  nonunique  solutions. 
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It  is  of  interest  to  note  that  the  foregoing  analysis  is 
based  on  the  boundary  conditions  of  Eguations  (6)  which  are  in 
similarity  form.  Hence  a  nonunique  solution  can  appear  for  a 
range  of  thickness  and  camber  parameters  provided  the  Mach  number 
is  changed  so  that  the  similarity  boundary  conditions  are 
unchanged.  This  is  consistent  with  the  numerical  evidence  of 
Salas  et  al .  (Ref.  2). 

2.7  Non-Conservative  Algorithms 

A  non-conservative  algorithm,  such  as  that  of  Murman  and 
Cole  (Ref.  7)  adds  source  terms  at  the  shock  waves.  Hence  the 
algorithm  is  solving  conservatively  a  differential  equation  of 

(l-<j>  )  <j>  +  4>  =  [o.H(x  -  x  )]  (60) 

Yx  Yxx  Yyy  l  x 

where  o.(y)  is  a  source  term  at  the  shock  location  x_  and 
1  &  l 

H  (  )  is  the  step  function.  The  strength  of  the  source  term  is 

not  known  explicitly  in  non-conservative  algorithms. 

If  the  asymmetric  part  of  Equation  (60)  is  taken  then 

A  +  A  $  =  (Auu)  +Mo.H(x-x  )]  (61) 

xx  yy  x  i  x 

If  Equation  (61)  is  added  to  Equation  (44)  the  following  result 
is  obtained. 

(A<J>+G)  +  (a<(>+G)  =  (Auu)  +  {g  +  A[o.H(x-x  ]}  (62) 

XX  I  i  **  ^  ^  ** 

If  a  ficticious  positive  component  of  lift  starts  to  appear 
in  a  solution,  g  will  change  from  zero  to  a  predominantly 
positive  quantity,  since  g  is  equivalent  to  (Auu).  This 
introduction  of  positive  lift  will  change  the  location  of  the 
shock  waves  such  that  an  upper  surface  shock  moves  aft  while  a 
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If  there  is  a  real  lift  due  to  the  boundary  condition  of  Equation 
(6)  then  Equation  (24)  gets  replaced  by 


/s,^tyl£+£*  -  L  1571 

where  I_  is  qiven  by  Equation  (28)  and  f*  is  a  function  that 
CL 

satisf ies 


L  /  ff*  -  I  ]d5dn  =  0 
bl  CL 


(58) 


and 


* 


The  analysis  described  above  is  unchanged  if  f  is  replaced 

•k  k 

by  f  +  f  since  f  always  occurs  as  part  of  the  term  f  +  f  ?  in 
this  case  f  +  f*  must  satisfy  Equation  (57)  rather  than  Equation 
(24).  Also,  since  f*  alone  will  satisfy  Equations  (56)  and  (57) 
there  must  be  at  least  two  consistent  solutions  for  f  +  f*  for 
multiple  solutions  to  exist. 


Equation  (56)  can  be  written  as 


l  +  I,  +  lt  +  g)  2  (q  +  lL  +  If+  Q)  =  ULS  +  ULA  4 n  ^S,-^Kx_[I£  + 


where 


+g+If)  d^dn 


+g 


(  59) 


LA 


-  /i 


-f  /(l-u  )  y 
o _ o 

2  2 

(x-c )  +  y 


dC 


Since  the  addition  of  the  uLA  term  has  the  same  effect  as 
introducing  camber  into  the  airfoil  a  nonunique  solution 
generally  will  have  a  different  shock  location  than  a  real 
solution . 


-14- 


it  also  follows  from  Equation  (53)  that 

AU  =  g/_  =  g  +  if 

From  Equation  (34) 


(53) 


(54) 


AU  = 


fo 


°  u-1 

o 

and  if  the  limit  as  y  +  0  is  taken  of  Equation  (31)  and  if 

Equations  (53)  and  (54)  are  used  then 

-2  2 

“o  -  5s  ULS  +  I  -  h  )is  U  /k«1(^T7)2  +  l9+if  >2)<3Wn 

O  O  I  t 

(  55) 


where  uLg  is  the  value  of  Uiq  at  y  =  0. 


LS 


Since  f  is  assumed  known  and  ULS„  is  known,  Equation  (55)  is  an 
equation  for  uq.  Hence  g  can  be  found. 


rising  Equations  (58)  and  (59)  to  eliminate  u  and  Au  in 
Eauation  (33)  it  follows  that  a  nonunique  solution  only  exists  if 
there  is  a  solution  f,  to  the  equation 


q+ 1 . 


1  ( — 9_ 

2  '  g  + 1 , 


)2  = 


LS 


(g+if ) 


L  I_  / 

4  77 


/Kx£  [  (g  +  I f ) 


2  +  (g  +  IfP]d£dn 


(56) 


This  solution  must  also  satisfy  Equation  (24).  For  a  real 
nonlifting  case  a  nonunique  solution  will  appear  if  a  nonzero 
value  of  f  satisfies  both  Equations  (24)  and  (56). 
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(49)  are  similar  in  form  to  Equations  (1)  and  (6)  if  $  is 
identified  with  $  +  G.  .  The  numerical  algorithm  will  solve 
Equation  (48)  in  an  identical  manner  to  the  real  solution. 

If  there  is  real  lift,  that  is  a  value  that  is  calculated 
with  no  eigensolutions ,  then  the  lifting  solution  of  Equation 
(38)  can  be  added  to  Equation  (48)  to  give 

( <(>  +  Al^+G^xx  +  ( ♦  “  A4l+G)yy  ~  \  +  Au+cj)2]^ 

u 

where  denotes  a  real  lifting  component.  In  this  case 

the  right  hand  side  of  Equation  (41)  is  identified  with 
[ Au  +  g/-] . 

The  boundary  condition  is 

1  (x,±0)  +  A  <}>*  ( x  ,  ±0  )  +  G  (x,±0)  =  Yg  ( x  ,±0  )  -  A  (51) 

Again  a  ficticious  lifting  component  G(x,y)  is  added  to  the 
equation  without  a  change  in  the  boundary  conditions.  If  G  is 
present  the  numerical  algorithm  will  solve  Equation  (50)  in  an 
identical  manner  as  the  real  solution.  It  is  not  possible  to 
decouple  the  real  and  ficticious  components. 

In  the  preceeding  analysis  u  has  been  assumed  known.  For  a 
solution  to  exist  it  also  must  be  a  solution  of  Equation  (31). 

Equation  (37)  can  be  written  as 

<^)  5  -  -  IT  lsj*tx  fd5dn  ’  Jf  (52> 

where  g  is  given  by  Equation  (46).  It  follows  from  Equation  (52) 
that 


(50) 

Au  on 
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Now  Equation  (37)  can  be  modified  to  give  the  integral  equation 


a  _  n  .  i_  r  ° 

—  ®  2  TT  '<*>  .  .  2  2 

u  (x-c)  +  y 


fo/(l-u  )y 

- - 75 - 5**  ~  77  U  /Kfvgd?dn 


4 it  JS^‘'Cx- 


which  can  be  written  in  differential  form  as 


G  +  G  =  g  (44) 

xx  yy  yx  '  ' 

where 

G(x,y)  =  -  g  U.y)  dC  (45) 

u ( s,y ) 

and 

q(C,y)  =  f(5,y) - (46) 

1-u ( 5»o) 

with  the  boundary  condition 

G  (x,±0)  =  n  (47) 

In  the  far  field  G  behaves  like  a  point  vortex. 

If  G(x,y)  is  identified  with  a  lifting  term  A<j>  then 
Equations  (41)  and  (44)  can  be  added  to  give 

<♦  +  G)xx+  +  G)yy  =  1  +  ai2x  (48) 

u 

ly(x,±0)  +  G  (x,±0)  =  -  [Yg(x,+0)-Yg(x,-0) ]  (49) 

where  Au  on  the  right  hand  side  of  Equation  (41)  has  been 
identified  as  g/u.  Thus  a  ficticious  value,  G,  can  be  added  to 
a  purely  symmetric  problem,  denoted  by  to  give  lift;  the 

boundary  conditions  are  not  affected.  This  is  the  mechanism  of 
the  appearance  of  the  nonunique  solutions.  Equations  (48)  and 
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(32).  Hence  for  a  nonunique  solution  to  exist  Au,  as  defined  by 
Equation  (35)  must  be  compatible  with  Equation  (32). 


Substitution  of  Equation  (33)  into  Equation  (32)  gives 

“u-4"o  ■  f  -  h  !SllKixc  d5d"  ,36: 

or,  using  Equations  (34),  (35) 

I  [f  .  fo  ii£l]  .  £  -  ;  /K  £  awn  (37: 

u  (1-uo)  1 

The  basic  differential  equation,  (1),  can  be  decoupled  into 
symmetric  and  asymmetric  parts.  The  asymmetric  part  is 


A  <b  +  A  di  =  (Auu) 
Yxx  Yyy  x 


where 


A<j>(x  ,y )  -  2  t<t>(x,y)  =  4>  ( x ,  —  y ) } 


the  boundary  condition  is 


A  <jy  ( x  ,  o )  =  -A  +  ~  [  Yg  (  x  ,  +0  )  +  (  x  ,  -0  )  ] 


The  symmetric  part  is 


—  1  -?  2 
j,  +  4,  =4  (u  +  Au  ] 

Kxx  Yxx  2  1  x 


with 


,(x,o)  =  \  [  Y^  (x  ,  +0  )  -  Y^  ( x ,  -0  )  ] 
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Equation  (10)  can  be  manipulated  to  give  the  following  symmetric 
and  asymmetric  parts. 

u(x,y)  -  ^  (u2(x,y)  +  Au2(x,y)]  =  UL  “  /s  /K?xtu2(C»n)  (31) 

+au  2  ( £, n) ] d?dn 

Au(x,y)  -  Au(x,y)u(x,y)  =  ~  ^  /s  [u2 ( n)  -  u2(C»n)  (32) 

-2 Au ( ? ,o ) ] dCdn 

These  equations  will  he  used  in  the  following  sections. 

2.6  Analysis  of  Multiple  Solutions 
If  in  Equation  (24),  f(C,n)  is  known  then  Equation  (23) 

gives 

~  [u2U,n)  =  u2U,n)-2AuU,0)]  =  Au(  C»n)u(  C,n)  “  au(C,0)  -fU,n) 

(33) 

If  iJ(5,ti)  is  known  then  this  is  an  equation  for  Au(£,n)  in  the 
flow  field  and  gives 


Au  =  (34) 

°  u  -1 
o 

whe  re 

f  =  f ( 5,0)  etc. ; 

also 

f+  Au  ,  r 

Au  =  -  o  m  1  ff - ^£_]  (35) 

u  u  1-u 

o 

In  the  transonic  solution  Equation  (16)  only  gives  the  value  of 
Au(£,o)  or  Auq,  the  value  of  Au(x,y)  being  found  from  Equation 


-9- 


where  for  0  <  x  <  1 


(x)  = 


1  -  x,  1/2 


/i 


-A 


T  1 


i  [YsU) 


Y*U)] 


(x  -  £) 


1  - 


1/2 


<U 


(28) 


and 


I  (x)  =0;  x<0,  x  >  1 ? 

CL 

u( C»n)  is  a  transonic  solution  subject  to  the  arbitrary  boundary 
conditions  in  Equation  (6).  Since  the  boundary  condition  is 
arbitrary,  and  therefore  an  infinite  number  of  boundary 
conditions  can  be  applied,  it  follows  that  an  infinite  number  of 
solutions  u(£,n)  exist. 

Equation  (28)  is  identical  in  form  to  Equation  (24)  with 
f  (  £, n)  qiven  by. 

fU,n)  =  u2U,n)  -  u2U,-n)  -  2au(£)  -21  (5)  /  0  (29) 

L 

Since  there  are  an  infinite  number  of  transonic  solutions  it 
follows  that  there  are  an  infinite  number  of  eigensolut ions 
f(£,n).  A  justification  of  a  nonzero  f (  £ , n )  is  given  in  the 
Appendix. 

2.5  Symmetric  and  Antisymmetric  Integral  Equations 
Let  u  (£,n)  and  Au(£,n)  be  defined  by 

u  (  £,  n)  =  ~  [u(£,n)  +u(£,-n)]  (30a) 

Au(£,n)  =  ^  (u(£,n)  +  u(£,-n)l  (30b) 
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(23) 


Jg/k^y  (u2  (  5,  n)  =  u2U,n)  -  2Au(  £)]d£dn  =  0 

where  covers  half  of  the  domain  S.  Equation  (23)  can  be 
written  in  the  more  compact  form. 

/g/lk^yf (C,n)dCdn  =  0  (24) 

where 

fU,n)  =  [u2U,n)  =  u2U,-n)  =  2au(£)] 

One  solution  of  this  integral  equation  is 

fU,n)  =  0  (25) 

which  is  the  symmetric  solution  for  the  airfoil  problem.  A 
nonsymmetric  solution  can  be  obtained  if  there  is  one  or  more 
functions  f(£,n)  =  /  0  that  satisfy  Equation  (24).  If  there  are 
such  solutions  then  a  multiple  solution  in  transonic  flow  can 
exist  if  the  equation 

\  (u2(£,n)  =  u2 ( £ , -n)  -  2Au ( £ ) ]  =  f(£,n)  (26) 

has  a  real  solution  for  u ( £ , n ) . 

2.4  Existence  and  Nature  of  Eigensolutions 

It  is  desirable  to  determine  the  existence  and  the  nature  of 
the  eigensolutions,  f(£,n)  of  Equation  (24).  Consider  Equations 
(16)  and  (18).  Using  the  inversion  procedure  that  leads  to 
Equation  (20),  Equations  (16),  (18)  can  be  written  in  the  form 

/  /K  [(u2(£,n)-u2(£,-n)  -2 Au ( £ )  -21  (£)]d£n  =  0  (27) 

al  c  L 
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This  integral  equation  is  valid  for  shock  waves  normal  to 
the  freestream;  if  the  shocks  are  not  normal  to  the  freestream  a 
modified  integral  equation  is  used  (Ref.  6).  Since  the 
formulation  changes  are  negligible  the  above  set  of  equations 
will  be  used  in  the  subsequent  discussions  for  clarity. 


In  the  solution  of  the  integral  equations  the  circulation  is 
given  by  Equation  (16)  although  it  should  be  emphasized  that 
Equations  (10)  and  (16)  are  not  independent.  The  solution  of 
Equation  (16)  is  found  by  inverting  the  integral  equation  to  give 


Au  ( x )  - 


Au  ( x ) 


=  u 


LA 


U.0)  -i(i^)1/2  /i 


1  o  x  -  5 


(r-H)1/2ac 

(20) 


where  uLA  is  the  antisymmetric  solution  of  Equation  (1)  without 
the  nonlinear  terms  and  is  given  by 


ULA{X'0) 


1/2  r  1 
1  O 


-A  +  \  [Yg(C,+0)  +Y^U,-0)] 
______ 


1/2 


dfi 


(21) 


The  inversion  procedure  invokes  the  Kutta  condition. 

2.3  Application  of  the  Integral  Equations  to  Multiple 
Solutions 

Consider  for  the  moment  the  case  of  a  symmetric  airfoil  at 
zero  angle  of  attack  in  this  case 

A  =  0  (22) 

Y  (x ,+0 )  =  -Y  (x , -0 ) 

S  5 

and  the  left  hand  side  of  Equation  (16)  is  zero.  A  manipulation 
of  the  integrals  in  Equation  (16)  leads  to  the  equation 


whe  re 


AuU)  =  [u  (  C  ,  +0 )  =  uU,-0)]/2  (13) 

The  Kernel  function  K(x,5;y»ri)  is  given  by 

K(x,£ry,n)  =  |  In  t(x  -  £)2  +  (y  -  n)2]  (14) 

The  integral  IT  is  given  by 

IT(u)  =  -  JSJ KCx(x,C?y,n)u2(C,n)dS  (15) 

The  domain  S  is  shown  in  Figure  2. 


If  y  =  ±  0,  Equation  (10)  gives  (see  Ref.  5)  only  the 
symmetric  part  of  the  solution  and  the  antisymmetric  part  is 
given  by 


-A  +  ~[Y£(x,±0) 


+  Y'(x,-0)}=  - 
s 


IAU ( S.?-x )_/!].  +  ic(X) 

(16) 


where 


AU2(U  =  (u2  (  ?  ,  +0  )  -  u2U,-0)]/2 


(17) 


and 


Ic(x)  =  -  JSJ  kKy (x,g;0,o>  [u2 ( C, n)  -  G2(5)]dS  (18) 


where 


*2 

uU) 


u2 ( C , +0 )  ,  n  >  0 
u2U,-0) ,  n  <  0 


( 19) 
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is  satisfied  and  that  the  Kutta  condition  of  zero  velocity  jump 
at  the  trailing  edge  and  on  the  wake  is  satisfied.  In  Equation 
(6) 


Yg ( x , ±0 )  =  k/B3Yg  ( x  ,  ±0 )  (7) 

where  y  =  Y  (x,±0)  denotes  the  geometry  of  the  airfoil  on  its 
upper  and  lower  surfaces,  respectively.  A  is  given  by 

A  =  k/83  a  (  8  ) 

where  a  is  the  angle  of  attack,  although  it  should  be  noted  that 
in  the  formulation  used  in  TSFOIL 

A  =  k/B3M  -1/4a  (9) 

00 

The  basic  idea  of  the  integral  equation  method  is  to  use  Green's 
theorem  to  write  the  differential  equation,  Eauation  (1),  and  its 
associated  boundary  conditions  in  integral  form.  A  detailed 
description  of  the  method  is  given  in  Reference  5. 

For  y  ^  0  the  integral  equation  is  given  by 

U  "  u2/2  =  ULS  +  ULA  +  1T{U)  (10) 


whe  re 


ULS(x,y)  =  hi  ^0[YsU,+0)  ~YsU'_0)]  Kx(x,C;y,0)dc  (11) 

and  is  the  solution  of  Eguation  (1)  without  the  nonlinear  terms, 
uLA(x,y)  =  I  /J  Au  (OK  (x,e ;y,0)d?  (12) 
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stability  terms  used  in  algorithms  tend  to  initiate  the 
appearance  of  e iqensolut ions  thus  making  the  computations  of  a 
real  solution  different. 

Since  the  behavior  of  the  nonunique  solutions  is  identical 
to  that  of  ^he  real  solution  the  question  arises  as  to  whether 
these  solutions  are  physically  realizable.  An  analysis  of  the 
Navier-Stokes  or,  possibly,  the  Euler  equations  is  necessary  to 
determine  this. 
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APPENDIX 


Further  Comments  On  The  Nature  of  Eigensolut ions 


It  is  possible  that  the  only  solution  to  Equation  (27)  is 


f(5,n)  =  0 


In  this  case 


4u  u  -  iu  -  I  =0 
°  CL 


CL 

Au  =  - 

°  u  -1 

o 


AU  =  -  I 


CL  u  -1 


If  Equation  (A2)  is  substituted  into  Equation  (32)  it  follows 
that 

-  i  i  :clU)  y 

Au  -  Au  u  =  —  /  - x - X  d£  ( A4  : 

11  J°  (x-o2+  y2 

The  right-hand  side  of  equation  (Ab)  can  be  recognized  as 
Au^,  the  linear  value  of  Au  in  the  flow  field.  Hence  Equation 
(A5)  becomes 


Au  ( 1  —  u )  =  Au 


It  can  be  seen  that  as  u  1  Au  becomes  infinite  which  will 
qive  an  expansion  shock.  Consequently,  Equation  (Al)  cannot  be 
the  correct  solution  of  Equation  (32);  a  nonzero  value  of 
f(f,,n)  must  exist. 
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APPENDIX  3 


Velocity  Induced  by  a  Moving  Blade 

The  object  is  to  estimate  the  induced  upwash  at  some  station 
y  due  to  the  oscillatory  motion  of  a  blade  at  y  =  0.  It  is 
assumed  that  y  is  sufficiently  far  from  the  oscillating  blade 
that  linear  theory  can  apply. 


From  classic  subsonic  theory  the  velocity  potential  in  the 
far  field  is  given  by  (see  Reference  15) 


Mx»y)  = 


(1) 


where 


r 


(x2  + 


?2)1/2 


(2) 


and 


k  =  ~~  ,  y  =  d-M2)1/2y 

1-MZ 


(3) 


Far  from  the  moving  blade  x  can  be  taken  as  zero  and,  on 
differentiation  with  respect  to  y. 


v(o ,y ) 


A  eilt//4,-l  iK  ,  iKy 

d-M2)!/2  k1/2  V3/2"  yixi 


(4) 


as  y  +  <*>,  v(o,y)  can  be  approximated  by 


v(o ,y ) 


-  iA 


( 1 — M  2 ) 1/2 


.  e'iKy 
“  A  =^72 


j<y2 

-1/2 

y 


eiKy  ein/4 


where  A  is  a  complex  constant. 


(5) 
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APPENDIX  4 


by  Alfred  Ayoub 


Convergence  of  the  Series 


ma 


n=l  n 


172 


l 


ina 


The  convergence  of  the  series  nl/2  Is  ^ue  Prlmarily  to 


ina 


as  the  series 


the  cancellation  effects  provided  by  the  tejijj  e 

y  1_  y  £ _ a 

L .  1/2  is  divergent  and  therefore  L,  1/2  is  not  absolutely 

n=l  n  n= 1  n 

convergent . 

For  values  of  c.  for  which  an  integer  p  exists  such  that 
“  ma  59 

y  £ _ 

pa  =  it,  the  series  1/2  can  be  reduced  to  an  alternating 

n=l  n 

series  by  grouping  all  consecutive  terms  of  the  same  sign;  the 
series  that  follows  has  monoton ically  decreasing  terms  and  is 
therefore  convergent. 


While  a  somewhat  similar  approach  can  be  followed  for  an 
arbitrary  a,  the  procedure  tends  to  become  auite  complex.  The 
following  proof  of  convergence  on  the  other  hand  is  more  direct 
and  exploits  equally  well  the  alternating  character  of  eina. 


ina 


n  =  l  n 


1/2 


l  a  .  e 
e  +  — 


i  2a 


ma 


1/2 


n 


1/2 


, ,  1  ,  i  a  ,  .  1  1  .  .  i  a ,  i2a. 

=  (1~  2 1/2  6  2 1/2  ”  3 1/2  (e  ] 


+  .  .  .  + 
r  1  1 


n1/2  (n+l)1/2 


,  ,  la  ,  ,  .ina,  , 

]  (e  +  . . .  +  e  )+... 
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<*>  in  a  00  . 

°r  E.  1/2  =  ^  [T/2 

n=l  n  n=l  n 


( n+1 ) 


i  ,  l  a  , 

J  (e  +  . 


.  in  a. 
.  +  e  ) 


The  series  on  the  right  hand  side  of  the  above  identity  is 
absolutely  convergent  since: 


n1/2  (n+1)1/2 


/  lot  . 

(e  +  • 


.  +  einot)  < 


1  /2 
2nJ/  z 


00 

r  1 

l  — is  convergent  (by  Cauchy's  integral  test] 
>«1  2nJ/^ 


.  ina  .  2 

and  I  ize  j  =  sin  n«/2_ 

1-e  a  sin2a/2 

<  - 1 -  finite  and  independent  of  M 

sin  a/2  for  a  *  o,  2tt  ,  .  .  . 

It  follows  that  the  series 


co  ,  , 

v  r  1  1  ,  , _ia  ,  ,  .ina, 

[  1/2  “  ,  .  >  l/2]  (S  ®  ) 

n=l  n  ( n+1 ) 


is  converge 


“  „  ina 

r  e 


nt  and  hence  £  T7>  1 

n=l 


s  convergent, 


It  is  interesting  to  note  that  the  same  proof  above  applies 

CO  CO 

to  all  series  of  the  form  J  a  b  where  I  b  is  bounded  but  not 

%  n  n  i  n 

n=l  n=l 

necessarily  convergent  and  the  terms  an  are  all  positive  or  all 
negative  and  an  tends  to  zero  monotonically  as  where  r  is  any 
positive  exponent. 
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